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ABSTRACT 


Consider a box that is filled with an ideal gas and that is aligned along 
Cartesian coordinates (x, y, z) having until length in the “y” direction and 
unspecified length in the “x” and “z” directions. Heat is applied uniformly 
over the “hot” end of the box (y = 1) and is removed uniformly over the 
“cold” end (y = 0) at a constant rate such that the ends of the box are main- 
tained at temperatures T 0 at y - 0 and T ] at_y = 1. Let U, V, and W denote the 
respective velocity components of a molecule inside the box selected at some 
random time and at some location (x, y, z). If T 0 = T v then U, V, and W are 
mutually independent and Gaussian, each with mean zero and variance RT 0 , 
where R is the gas constant. When T 0 * T, the velocity components are not 
independent and are not Gaussian. Our objective is to characterize the joint 
distribution of the velocity components U, V, and IT as a function ofy, and, in 
particular, to characterize the distribution of V given >>. It is hoped that this 
research will lead to an increased physical understanding of the nature of 
turbulence. 
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THE SIMULATION 



Consider a box that is filled with an ideal gas and that is aligned along Cartesian coor- 
dinates ( x , y, z) having until length in the “y” direction and unspecified length in the “x” and 
“z” directions. Heat is applied uniformly over the “hot” end of the box (y = 1) and is re- 
moved uniformly over the “cold” end (y = 0) at a con- 
stant rate such that the ends of the box are maintained 
at temperatures T Q at y = 0 and T, at y = 1 . Let U, V, 
and W denote the respective velocity components of a 
molecule inside the box selected at some random time 
and at some location (x, y, z). (See Figure 1 .) Our con- 
cern in this project is to characterize these velocity 
components. In this regard, we will first consider the 
problem of simulating the molecular motion inside the Figure 1 

box so as to obtain simulated velocities. This simulation will be based on the direct simula- 
tion Monte Carlo code presented in [Bird, 1994], 

For the simulation, we will consider (as indicated by Figure 2) a thin slice of the box 
that effectively ignores the z spatial component. The z component is not entirely neglected 
however since the molecule velocities are considered three dimensional even though the 
molecule positions are taken to be only two dimensional. This two dimensional slice is di- 
vided into cells and subcells as indicated by Figure 3. As Bird describes the problem in [Bird, 
1994, p. 215], subcells were introduced because it “was feared that coincidental collisions 
between molecules at opposite sides of the cells, coupled with the immediate migration of 
these molecules to the next cell, might cause false disturbances to propagate with speed 

equal to the ratio of the cell 
size to the time step.” 

In order to explain the 
operation of the simulation 
program, we will consider the 
input file dsmc2 .dat that 
contains the data used to pro- 
duce a particular simulation 
run. Any line in the data file 
beginning with *#’ is treated 
as a comment and is ignored 
by the program. The line numbers mentioned below refer to the copy of this file that is 
reproduced at the end of this report. 
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The first item in the file (lines 6-9) determines whether or not an output file dsmc2 . out 
is produced. This file corresponds to the only output provided by the original program 
dsmc2 .for from [Bird, 1994], 

When the program begins, it first asks whether a new calculation is desired or an old 
calculation should be restarted. This latter option is only possible when line 1 9 of the data file 
has been set to 0 in an earlier run in order to produce a restart file. The restart file dsitic2 .res 


contains a snapshot of the current state 
of the simulation and can be used to re- 
start the program at that point should 
any sort of interruption occur. (Note that 
the size of this restart file can be as large 
as 20 Meg.) For long simulations, it is 
suggested that the values in lines 19-20 
be chosen so that a restart file is pro- 
duced approximately once each hour. 
Producing the file more frequently may 
greatly reduce the speed of the simula- 



Figure 3 


tion. 

A second output option is the production of alog file dsmc2 . log that contains veloc- 
ity and temperature information for cells along a particular y location. This file is produced 
when line 26 of the data file is set to zero. The value for^ is specified by line 1 50 of the input 
file. In addition, this log file contains information regarding the number of molecules per cell 
that may be useful is choosing a proper value of CWRY, which is considered below. A sample 
of this output file is given near the end of this report. Note that the averages given by this 
option are not cumulative but are instead based only on the number of molecules occupying 
cells along the given .y location at the time the file is updated. 

The next output option determines whether or not a chart of velocity moments is pro- 
duced during the simulation. When this option is selected by setting line 3 1 of the data file to 
zero, a chart as given later in this report is periodically written to the screen and to the file 
dsmc2 . cht. The statistics in this chart are cumulative over the run of the simulation. The 
third column indicates the total number of molecules that have been used to obtain the sample 
moments found in the final three columns. The row numbers that are used to produce the 
chart are determined by lines 163-182 of the data file. Note that it is possible to produce this 
data for every row in the box simply by setting line 167 in the data file equal to the total 


number of rows in the box. 

A third output option makes it possible to track the position of a molecule initially 
located along some specified row. A file called dsmc2 .mol containing the (x, y) coordi- 
nates of a molecule in the specified row is produced when line 36 of the data file is set to 
zero. The initial molecule row number is specified by line 1 54 of the data file. Since the 
molecules generally move only slightly between successive time steps, the step size provided 
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by line 161 makes it possible to more easily obtain a picture of the molecule’s motion over a 
longer period of time.. A graph that shows how the data file dsmc2 .mol can be used to 
study the molecular motion produced by the simulation is given below. 

The cell widths in they' direction are not uniform, but instead are chosen to be smaller 
at the cold end of the box than at the hot end of the box. The ratio CWRY of the cell width 
at the hot end of the box to the cell width at the cold end is given in line 40 of the data file. If 

S = CWRF" 99 , then the . . 

width of the cell in Row / is 

given by WS*~ l ( 1 -S)/( 1 -S 200 ) ' / 

where W is the total width „. MM . \. 

of the box in the y direction. — ± = — 

Note that here we have as- " "" 

sumed that the total number o \ 

of rows in the y direction is 

200. This value can be >•«« • — — ' j / 

changed if the code is c — 

recompiled. Generally, the 

ratio CWRY should be the 1 

C.0CS-! D. 05572 0.00514 0.CO37C 0.005TB 0.005B 0.00502 0.005B4 0.0Q5B6 0.50508 0.0059 0.9059: 

square root of the ratio of 

the temperatures at the two ends of the boxes. A graph showing the number of molecules in 
each row for a particular simulation is given near the end of this report. 

The initial temperature of the gas is determined by FTMP, which we have chosen to be 
the geometric mean of the two surface temperatures. This value is given in line 44 of the data 
file. The initial number density FND is in molecules per square meter. This value is specified 
in line 48 of the data file. 

In the DSMC method, a single simulated molecule represents FNUM actual molecules. 
This value is given in line 52 of the data file and must be chosen so that the total number of 
» x simulated molecules is not greater 

CBf2 j ounoary than the maximum allowed num- 

ber of molecules, which in this case 

. . is 40000. (This value can be 

Boundary 3 Boundary 4 . 

3 changed if the code is recompiled.) 

y The actual number of simulated 

CB(1] 1 » molecules used by the program is 

Boundary 1 given by ( FND/FNUM ) multiplied 

cb[3] cb[4] by the area of the box. 

Figure 4 The value DTM in line 56 of 

the data file is the time step of the 
simulation, and it should be chosen so that the molecular motion and the collisions are un- 
coupled. In particular, it should be much less than the local mean collision time. The DSMC 


Boundary 3 


Boundary 4 


Boundary 1 


Figure 4 
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method becomes increasingly more accurate as the cell size and the time step tend to zero. 
(See [Bird, 1994, p. 215].) 

The values CB[i ] for / = 1, 2, 3, and 4, are specified in lines 69-81 of the data file and 




Surface 1 lies along 
lower boundary of 
Row 1 



Row 200 



Surface 2 lies along 
lower boundary of 
Row 201 


Figure 5 


determine the dimen- 
sions of the box as indi- 
cated by Figure 4. Two 
planar surfaces (Surface 
1 and Surface 2) at pos- 
sibly different tempera- 
tures may be placed any- 
where within the box. 
The reflection at these 
surfaces is diffuse, which 
means that the return ve- 
locity of a colliding mol- 
ecule is determined by 
the temperature of the 


wall instead of by the incident velocity. At the walls (or boundaries) of the box, the reflection 
is specular, which means that the appropriate velocity component before the collision is 
simply reversed to obtain the velocity component after the collision. 

Surface 1 is specified by the parameters in lines 92-110 of the data file. The value 
ISURF[l] specifies the direction of a line normal to the surface and pointing toward the gas. 
For specular reflection, it is this velocity component that is reversed upon reflection. (The 
velocity components parallel to the surface are not changed.) For example, in Figure 5, such 
a normal is in the positive.)' direction, which corresponds to a value of 1 for ISURF[ 1]. The 
negative .y direction corresponds to a value of 2, the positive x direction corresponds to a 
value of 3, and the negative x direction corresponds to a value of 4. The values of LIM[l][i] 
specify the location of Surface 1 as indicated by Figure 5 and the descriptions in lines 96-106 
of the data file. The value TSURF[\] is the temperature of Surface 1. Surface 2 is similarly 
specified by the parameters in lines 112-130. 

The values in lines 132-146 determine the rate of output and the duration of the simu- 
lation. Their use is illustrated by the flow chart near the end of this report. Note that the 
sampling does not begin until NPS output cycles have occurred. The moment chart is not 
produced until this same point is reached. 


FUTURE DIRECTIONS 


The ultimate goal of this work is to characterize the conditional distribution of the.y 
component of the velocity given a particular location^. Since this problem appears to be 
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analytically intractable, any practical investigation of the distribution will involve data ob- 
tained either empirically or via a simulation such as the one that we have considered here. 

One immediate item for future consideration involves the effect of the random number 
generator on the simulation results. We, as did Bird, have used the algorithm by Knuth that 
is found in [Press etal., 1992, p. 283], In [Ferrenberg etal, 1992], it was shown that subtle 
numerical errors can occur in Monte Carlo simulations when many supposedly ‘good’ ran- 
dom number generators are used. Since our simulation requires a large number of random 
numbers, it would be helpful to repeat and compare these simulations with a variety of 
different random number algorithms. 

A second item for future consideration involves an investigation of the boundary condi- 
tions at the two surfaces. The current code assumes that the distribution of the velocity 
component that is normal to the surface has the form of a zero mean Gaussian distribution 
taken to the right of the origin and renormalized. This is the distribution that would corre- 
spond to molecules crossing the surface from some external stationary equilibrium gas. (See 
[Bird, 1 994, p. 259] .) An alternate model that provides a better fit with the simulated data is 
a distribution obtained by taking the absolute value of the difference of two independent 
gamma random variables. 
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10000 




Scatter plot ofy velocity component in meters/second for approximately 16000 molecules 
selected at random from the box after approximately 268 million collisions have occurred. 
The temperature at row 1 is 300 and at row 200 is 30000. 
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A plot of the root mean square value of they velocity component for each row in the box. 
Again, the values were obtained for a box that had experienced approximately 268 million 
simulated collisions. The temperature at row 1 is 300 and at row 200 is 30000. 
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Sample of dsmc2 . log 



Molecules 
Collisions 
Start Time 
End Time 
DTM 


Surf 1 Temp 
Surf 2 Temp 


37538 
145790602 
1 . 1 6e-06 
0.00055646 
2e-09 
1 . 3e+l 4 
10 

3000 

0 

300 

30000 


Row Y Midpoint 
20 0.00055534 

40 0.00127118 

50 0.00169623 

60 0.00217341 

80 0.00331058 

100 0.00474384 

120 0.00655029 

140 0.00882712 

160 0.01169679 

180 0.01531368 


Mol Count 
2874478 
2260561 
2158229 
2103919 
2086605 
2146617 
2246682 
2378198 
2527745 
2695610 


2nd Moment 
0.7194 
1.1495 
1.3561 
1.5588 
1.9812 
2.4251 
2.9223 
3.4810 
4.1265 


ird Moment 
-0.2035 
-0.3343 
-0.3945 
-0.4584 
-0.5631 
-0.7288 
-0.8890 
- 1.0843 
-1.3110 
-1.5899 


4th Moment 
1.7523 
4 . 3424 
6.0186 
7.9019 
12.6606 
18.9076 
27.4045 
38.6212 
53.6066 
73.5440 



9-12 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 


# 

# DATA FILE FOR DSMC2 

# Any line beginning with is ignored. 

# Data should be placed on a single line following each description. 

# 

# 1 : Do not produce an output file 

# 0: Produce output files 

# 

1 

# 

# 1 : Do not produce a restart file 

# 0: Produce a restart file at every Nth output file 

# The value of N should follow the flag value below. 

# For example, if the next two data lines were 

# 0 

# 5 

# then a restart file would be produced after each group of 5 outputs. 



0 

100 

# 

# 1 : Do not produce a log file 

# 0: Produce a log file whenever an output file is produced 

# See Y Value below 

# 

1 

# 

# 1 : Do not produce a moment chart 

# 0: Produce a moment chart whenever an output file is produced 

# 

0 

# 

# 1 : Do not produce a molecule track file 

# 0: Track a molecule in row specified below 

# 

0 

# 

# CWRY : Ratio of cell width at outer to that at inner y boundary 

# 

10 

H 

# FTMP: The stream temperature 

# 

3000 



# FND: The initial number density 



1 . 22e22 

# 

# FNUM: The number of real mols represented by each simulated mol 



0. 13el5 



# DTM: The time step over which mol motion and collisions are uncoupled 

I 

2e-9 



# Boundary 2 

# high x 

# I i 

# Boundary ! I Boundary 

#31 14 
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63 

64 

65 

66 

67 

68 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

79 

80 
81 
82 

83 

84 

85 

86 

87 

88 

89 

90 

91 

92 

93 

94 

95 

96 

97 

98 

99 
100 
101 
102 

103 

104 

105 
1 Co 

107 

108 

109 

110 
111 
112 

113 

114 

115 

116 

117 

118 

119 

120 
121 
122 

123 

124 


ft t 1 

ft low X 

ft low y Boundary 1 high y 

# 

ft CB[1]: The x coordinate of boundary 1 (low x) 

ft 

-0.01 

ft 

ft CB [ 2 ] : The x coordinate of boundary 2 (high x) 

* 

0.01 

ft 

ft CB[3]: The y coordinate of boundary 3 (low y) 

ft 

0 

ft 

ft CB [ 4 ] : The y coordinate of boundary 4 (high y) 

# 

0.02 

fi 

ft 

Col 20 

# I I 

ft Surface I > pos y neg y < 1 Surface 

ft 1 ! dir dir I 2 

ft ! I 

ft Col 1 

ft Row 1 Row 200 

# 

ft ISURFU1 : Direction of surface 1 normal is in the positive y direction 

ft 

1 

ft 

ft LIMS[1][1]: Surface 1 lies along lower boundary of this row 

ft 

1 

ft 

ft LIMS [ 1 ] [ 2 ] : Surface 1 begins at column 1 

ft 

1 

ft 

ft LIMS [ 1 ] [3] : and continues until column 20 

ft 

20 

ft 

ft TSURF[1]: Temperature of surface 1 

ft 

300 



ft ISURF [ 2 ] : Direction of surface 2 normal is in the negative y direction 

ft 

2 

ft 

# LIMS [2] [ 1 ] : Surface 2 lies along lower boundary of this row 

ft 

201 

ft 

ft LIMS [ 2 ] [ 2 ] : Surface 2 begins at column 1 

ft 

1 

ft 

ft LIMS [ 2 ] [ 3 ] : and continues until column 20 
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125 

126 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138 

139 

140 

141 

142 

143 
14 4 
145 
14 6 

147 

148 
14 9 

150 

151 

152 

153 

154 

155 

156 

157 

158 

159 

160 
161 
162 

163 

164 

165 

166 

167 

168 

169 

170 

171 

172 

173 

174 

175 

176 

177 

178 

179 

180 
181 
182 


If 

20 



ff TSURF[2}; Temperature of surface 2 

ff 

30000 



ft NIS: The number of DTM time steps between samplings 

ft 

2 



# NSP: The number of samples between prints 



10 

If 

If N PS : The number of prints until assumed start of steady flow 

If 

30 

ff 

If NPT: The total number of prints 



50000 



If Y value for which output will be collected in the LOG FILE 



0.OC584 

If 

If Row Number Containing Molecule to Track 



ICO 

If 

ff Step size for molecule position output to file 

# i.e. if next data line is 10 then position information will be 

# written to the file every 10 time steps 

# Use 1 to output a position after every time step 

If 

10 

ff 

ff Number of Rows for which Output is Desired. 

If If this value equals or exceeds the total number of rows then data 

# will be output for all rows. 



I 

| 

!* Row Numbers to Output. 

I # Required if previous data line is greater than zero and less than the 
| # total number of rows. 

| 

20 

I 40 

| 50 
60 
80 

I 100 
I 120 

|140 
160 
I 180 
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